doseresNMA antidep/drnma.model.standardFE.R

In htx-r/doseresNMA:

# FE delta's
drnma.standardFE <- function(){
  for(i in 1:ns){ # loop through ns trials
    delta[i,1] <- 0
    u[i] ~ dnorm(0,0.001) # baseline effect

    for (k in 1:na[i]){ # loop through na arms within a study
      r[i,k] ~ dbin(p[i,k], n[i,k])
      logit(p[i,k]) <- theta[i,k]
      theta[i,k] <- u[i] + delta[i,k]

      # arm deviance
      rhat[i,k] <- p[i,k] * n[i,k] # expected number of events
      resdev[i,k] <- 2 * (r[i,k] * (log(r[i,k]) - log(rhat[i,k])) + (n[i,k] - r[i,k]) * (log(n[i,k] - r[i,k]) - log(n[i,k] - rhat[i,k])))

    }

    # study deviance
    resstudydev[i] <- sum(resdev[i, 1:na[i]])

    for(k in 2:na[i]){ # treatment effects
      # FE delta's
      delta[i,k] <- DE[i,k]
      # dose-effect model
      DE[i,k] <- ((beta1[t[i,k]] * dose1[i,k]) + (beta2[t[i,k]] * dose2[i,k])) - ((beta1[t[i,1]] * dose1[i,1]) + (beta2[t[i,1]] * dose2[i,1]))
          }
  }

  #** PRIORS
  # beta's priors
  beta2[ref.index] <- 0
  beta1[ref.index] <- 0
  for (k in c(1:(ref.index - 1), (ref.index + 1):ndrugs)){ # Priors on relative treatment effects
    beta1[k]~dnorm(0, 0.001)
    beta2[k]~dnorm(0, 0.001)
  }

  # overall deviance
  totresdev <- sum(resstudydev[])

  #** drug absoulte effect
  # baseline effect
  for (m in 1:ns){
    rr[m,1] ~ dbinom(p0[m],nn[m,1])
    logit(p0[m]) <- OO[m]
    OO[m] ~ dnorm(Z, prec.z)
  }
  # baseline priors
  Z ~ dnorm(0, 0.001)
  prec.z <- pow(sigma.z,-2)
  sigma.z ~ dunif(0,10)
  p.placebo <- exp(Z)/(1+exp(Z))
  # combine baseline with OR -> drug absolute effect
  for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
    for( j in 1:nd.new){
      OR[j,k] <- exp(beta1[k]*new.dose[k,j]+  beta2[k]*f.new.dose[k,j])
      odds.drug[j,k] <- OR[j,k]*exp(Z)
      p.drug[j,k] <- odds.drug[j,k]/(1+odds.drug[j,k])
    }
  }
}